Meshless method for solid mechanics simulation, electronic device, and storage medium

ABSTRACT

A meshless method for solid mechanics simulation, which includes the following steps: first constructing a trial function and a test function by distributing points according to a geometric shape of a structure; then introducing numerical flux to correct the traditional Galerkin weak form; and finally substituting the trial and test functions into the weak form to construct a global stiffness matrix of the structure and a system of algebraic equations, and solving the system of algebraic equations to obtain displacement and stress of each point, thus completing simulation and analysis of the structure deformation and stress.

CROSS REFERENCE TO RELATED APPLICATIONS

The present application claims priority to the Chinese PatentApplication No. 201910629236.4 entitled “MESHLESS METHOD FOR SOLIDMECHANICS SIMULATION, ELECTRONIC DEVICE, AND STORAGE MEDIUM”, filed onJul. 12, 2019, and incorporates its disclosure herein by reference inits entirety.

TECHNICAL FIELD

The present disclosure relates to a mechanics simulation system, and inparticular, to a meshless method for solid mechanics simulation, anelectronic device, and a storage medium.

BACKGROUND

Structural stress analysis is very important and necessary in diverseengineering fields, such as aeronautics, astronautics, automobileengineering, etc. From design, manufacture to maintenance of products,structural stress analysis always plays a crucial role. Because of itssignificance, numerous researchers have been focusing on improving theaccuracy and efficiency of the process of structural stress analysis fordecades.

Finite element method (FEM) is widely used in structural stress analysisat present. This method resorts to element-based, local, polynomial, andcontinuous trial and test functions for simulation. Because the trialand test functions are all polynomials, integrals in the FEM are easy tocompute. Moreover, the symmetry and sparsity of a global stiffnessmatrix make the FEM suitable and efficient in large-scale structuralsimulations. However, the accuracy of the FEM highly depends on thequality of mesh. In order to obtain a satisfactory solution, a lot ofefforts are usually spent on meshing. Further, even if simulation isinitialized with a high-quality mesh structure, mesh distortion mayoccur in the FEM in the case of a large structural deformation,dramatically reducing the precision of the solution. Further, in orderto enable FEM to investigate the formation of cracks, and the ruptureand fragmentation of structures, it is necessary to repeatedly performre-meshing, element deletion or use a cohesion model. However, theseadditional means render the FEM time-consuming and laborious incalculating extreme problems, and the result is not accurate enough.

Meshless methods, which eliminate meshes partly or completely, have beendeveloped since the end of last century. The meshless method is mainlyclassified into two types: a meshless weak form method and a meshlessparticle method. Element-free Galerkin (EFG) and Meshless LocalPetrov-Galerkin (MLPG) are two typical meshless weak-form methods, whichmainly utilize Moving Least Squares (MLS) or Radial Basis Function (RBF)approximations to deduce point-based trial functions. The use of the MLSand RBF approximations makes it easy to achieve higher-order continuityof the functions. Besides, by replacing the mesh structure withindividual points, the meshless method can conveniently insert or removeadditional points. However, on the other hand, the trial functionsobtained by the MLS or RBF method are no longer polynomial functions,and are rather complex. Therefore, the integral calculations in both EFGand MLPG are extremely cumbersome, less accurate, and may even influencethe precision and stability of the whole method. Definitely, to reducethe computational cost and improve the accuracy of integration, somespecial, novel numerical integration methods, for example, a series ofnovel nodal integration methods, are necessarily adopted.

Smoothed Particle Hydrodynamics (SPH) method, as a meshless particlemethod, is simple in form and needs less computational cost. However,the SPH method is based on strong form rather than weak form. Thismethod has poor stability and is prone to problems such as tensileinstability. As another example, Peridynamics method is a relatively newmethod, but the theoretical basis thereof is different from thetraditional continuum mechanics. Classic constitutive models andengineering experiences accumulated by scientists and engineers fordecades are difficult to be directly used in the Peridynamics method.

SUMMARY

To overcome the shortcomings in the prior art, a first objective ofembodiments is to provide a meshless method for solid mechanicssimulation which can improve calculation efficiency and simulationprecision, and solve the problem that traditional numerical methods aredifficult to simulate some extreme problems such as rupture andfragmentation.

A second objective of the embodiments is to provide an electronicdevice, which can improve calculation efficiency and simulationprecision, and solve the problem that traditional numerical methods aredifficult to simulate some extreme problems such as rupture andfragmentation.

A third objective of the embodiments is to provide a computer-readablestorage medium, which can improve calculation efficiency and simulationprecision, and solve the problem that traditional numerical methods aredifficult simulate some extreme problems such as rupture andfragmentation.

The first objective is implemented through the following technicalsolution:

A meshless method for solid mechanics simulation includes the followingsteps:

a step of constructing a trial function and a test function, whereinobtaining expressions of the trial function and the test function withpoint displacements in each subdomain according to points distributed ina structure;

a step of developing a Galerkin weak form, wherein introducing numericalflux corrections to obtain a Galerkin weak form with the numerical fluxcorrections;

a step of generating and solving a system of algebraic equations,wherein constructing a global stiffness matrix of the structure and asystem of algebraic equations according to the expressions of the trialand test functions in each subdomain and the Galerkin weak form with thenumerical flux corrections, and solving the system of algebraicequations to obtain the displacement of each point; and

a step of calculating structural deformation and stress, whereincalculating the deformation and the stress of each position in thestructure according to the displacement of each point, to completesimulation and analysis of the structure.

In some embodiments, the method further includes a simulating step foraccording to a set sequence of loads in chronological order, performingthe step of constructing trial functions and test functions, the step ofdeveloping the Galerkin weak form, the step of generating and solving asystem of algebraic equations, and the step of calculating structuraldeformation and stress under each load successively, to obtain thedeformation and stress of each position in the structure; anddetermining, according to a certain criterion, whether there is a crackbetween neighboring subdomains in the structure, to simulate rupture andfragmentation of the structure under the effect of the sequence of loadsin chronological order.

In some embodiments, each subdomain contains one point.

In some embodiments, a process of constructing the trial function andthe test function is as follows:

step S11: first defining the trial function in each subdomain andrecording the trial function as u^(h), where the trial function is annth-order polynomial function, and n is a natural number greater than orequal to 1 and often equal to 1; and performing nth-order Taylorexpansion of the trial function at a point P₀, as shown in the followingformula (1):

$\begin{matrix}{{{u^{k}\left( {x,y} \right)} = {\begin{bmatrix}u_{x}^{h} \\u_{y}^{h}\end{bmatrix} = \begin{bmatrix}\begin{matrix}{u_{x}^{0} + {\frac{\partial u_{x}}{\partial x}{_{P_{0}}{\left( {x - x_{0}} \right) + \frac{\partial u_{x}}{\partial y}}}_{P_{0}}\left( {y - y_{0}} \right)} + \ldots +} \\{\frac{1}{n!}\frac{\partial^{n}u_{x}}{\partial y^{n}}\left( {y - y_{0}} \right)^{n}}\end{matrix} \\\begin{matrix}{u_{y}^{0} + {\frac{\partial u_{y}}{\partial x}{_{P_{0}}{\left( {x - x_{0}} \right) + \frac{\partial u_{y}}{\partial y}}}_{P_{0}}\left( {y - y_{0}} \right)} + \ldots +} \\{\frac{1}{n!}\frac{\partial^{n}u_{y}}{\partial y^{n}}\left( {y - y_{0}} \right)^{n}}\end{matrix}\end{bmatrix}}},} & (1)\end{matrix}$

where [u_(x) ⁰ u_(y) ⁰]^(T) is value of the displacement u^(h) at thepoint P₀, and

${{a = \begin{bmatrix}{\frac{\partial u_{x}}{\partial x}}_{P_{0}} & {\frac{\partial u_{x}}{\partial y}}_{P_{0}} & \ldots & \frac{\partial^{n}u_{x}}{\partial y^{n}} & {\frac{\partial u_{y}}{\partial x}}_{P_{0}} & {\frac{\partial u_{y}}{\partial y}}_{P_{0}} & \ldots & \frac{\partial^{n}u_{y}}{\partial y^{n}}\end{bmatrix}}}_{P_{0}}$is a displacement gradient from the first order to the nth order at thepoint P₀;

step S12: determining the support domain of each subdomain, and defininga weighted discrete L² norm J, where

$\begin{matrix}{{J = {\left( {{Aa} + u_{0} - u_{m}} \right)^{T}{W\left( {{Aa} + u_{0} - u_{m}} \right)}}};} & (2)\end{matrix}$ $\begin{matrix}{A = \begin{bmatrix}{x_{1} - x_{0}} & {y_{1} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{1} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\{0} & 0 & {\ldots} & 0 & {x_{1} - x_{0}} & {y_{1} - y_{0}} & {\ldots} & {\frac{1}{n!}\left( {y_{1} - y_{0}} \right)^{n}} \\{x_{2} - x_{0}} & {y_{2} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{2} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & {x_{2} - x_{0}} & {y_{2} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{2} - y_{0}} \right)^{n}} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\{x_{m} - x_{0}} & {y_{m} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{m} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & {x_{m} - x_{0}} & {y_{m} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{m} - y_{0}} \right)^{n}}\end{bmatrix}} & \end{matrix}$ $\begin{matrix}{{u_{0} = \left( \begin{bmatrix}{u_{x}^{0}} & u_{y}^{0} & u_{x}^{0} & u_{y}^{0} & \ldots & u_{x}^{0} & u_{y}^{0}\end{bmatrix}_{1 \times 2m} \right)^{T}},} & \end{matrix}$ $\begin{matrix}{{u_{m} = \begin{bmatrix}{u_{x}^{1}} & u_{y}^{1} & u_{x}^{2} & u_{y}^{2} & \ldots & u_{x}^{m} & u_{y}^{m}\end{bmatrix}^{T}},} & \end{matrix}$ $\begin{matrix}{{W = \begin{bmatrix}w_{x}^{1} & 0 & 0 & \ldots & 0 \\0 & w_{y}^{1} & \ddots & \ddots & \vdots \\0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & w_{x}^{m} & 0 \\0 & \ldots & 0 & 0 & w_{y}^{m}\end{bmatrix}},} & \end{matrix}$

(x_(i), y_(i)) is coordinates of a point P_(i), [u_(x) ^(i) u_(y)^(i)]^(T) is value of the displacement u^(h) at the point P_(i), and[w_(x) ^(i) w_(y) ^(i)]^(T) is value of a weight at the point P_(i),(i=1, 2, 3, . . . m).

step S13: expressing the displacement gradient a with the pointdisplacement according to the formula (2) and stationary conditions ofthe norm J, that is:

$\begin{matrix}{{a = {Cu}_{E}},} & (3)\end{matrix}$ $\begin{matrix}{{{{wherein}C} = {\left( {A^{T}{WA}} \right)^{- 1}A^{T}{W\begin{bmatrix}I_{1} & I_{2}\end{bmatrix}}}},} & \end{matrix}$ $\begin{matrix}{{u_{E} = \begin{bmatrix}u_{x}^{0} & u_{y}^{0} & u_{x}^{1} & u_{y}^{1} & \ldots & u_{x}^{m} & u_{y}^{m}\end{bmatrix}^{T}},} & \end{matrix}$ $\begin{matrix}{{I_{1} = \begin{bmatrix}{- 1} & 0 \\0 & {- 1} \\{- 1} & 0 \\0 & {- 1} \\ \vdots & \vdots \\{- 1} & 0 \\0 & {- 1}\end{bmatrix}_{2m \times 2}},{{{{and}I_{2}} = \begin{bmatrix}1 & 0 & \ldots & 0 \\0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\0 & \ldots & 0 & 1\end{bmatrix}_{2m \times 2m}};}} & \end{matrix}$

step S14: obtaining the expression of the trial function in eachsubdomain according to formulae (3) and (1); and

step S15: obtaining the expression of the test function by prescribingthe test function in the same form as that of the trial function.

In some embodiments, a definition of the support domain includes, but isnot limited to, any one of the following manners: defining the supportdomain of each point to contain all subdomains having points adjacent tothe corresponding point located therein, or defining the support domainof each point to contain all subdomains having points therein which arelocated in a circle of a particular size and with the correspondingpoint at the center; and if a crack exists between neighboringsubdomains, the point located on one side of the crack is not containedin the support domain of the point on the other side of the crack.

In some embodiments, the numerical flux corrections include, but are notlimited to, interior penalty numerical flux corrections.

In some embodiments, the step of generating and solving an system ofalgebraic equations further includes:

step S21: substituting the trial and test functions in each subdomaininto the Galerkin weak form corrected by numerical fluxes, to obtain apoint stiffness matrix and an internal-boundary stiffness matrix, wherethe point stiffness matrix is K_(E)=∫_(E)B^(T)DBdΩ, theinternal-boundary stiffness matrix is

${K_{h} = {{- {\int_{e}{\left\lbrack N^{T} \right\rbrack\left\{ {n_{e}{DB}} \right\} d\Gamma}}} - {\int_{e}{{\left\{ {B^{T}{Dn}_{e}^{T}} \right\}\lbrack N\rbrack}d\Gamma}} + {\frac{\eta}{h_{e}}{\int_{e}{{\left\lbrack N^{T} \right\rbrack\lbrack N\rbrack}d\Gamma}}}}},$and e is a common boundary between two neighboring subdomains;

step S22: obtaining a global stiffness matrix K by assembling all thepoint stiffness matrices and the internal-boundary stiffness matrices,to obtain the following system of algebraic equations in a matrix form:Kq=Q, where q is a point displacement vector, and Q is a load vector;and

step S23: solving the system of algebraic equations according to theglobal stiffness matrix and the load vector to obtain a pointdisplacement vector of the structure, thus finally obtaining adisplacement of each position in the structure.

The second objective is implemented through the following technicalsolution:

An electronic device includes a memory, a processor, and a computerprogram stored in the memory and executable in the processor, where theprocessor implements the steps in the above meshless method for solidmechanics simulation when executing the program.

The third objective is implemented through the following technicalsolution:

A computer-readable storage medium which stores a computer program isprovided. When executed by a processor, the computer program implementsthe steps in the above meshless method for solid mechanics simulation.

The present disclosure has the following advantageous effects comparedto the prior art:

The present disclosure forms a novel meshless method for solid mechanicssimulation by providing new discontinuous trial and test functions and acorresponding weak form, which can improve calculation efficiency andsimulation precision, and solve the problem that traditional numericalmethods are difficult to simulate some extreme problems such as ruptureand fragmentation.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a schematic diagram of a structure and points randomlydistributed therein provided by the present disclosure;

FIG. 2 is a schematic diagram of the structure and division ofsubdomains thereof provided by the present disclosure;

FIG. 3 is a schematic diagram of a support domain of point P₀ providedby the present disclosure;

FIG. 4 is a schematic diagram of a discontinuous shape function providedby the present disclosure;

FIG. 5 is a schematic diagram of a discontinuous trial function providedby the present disclosure;

FIG. 6 is a flowchart of a meshless method for solid mechanicssimulation provided by the present disclosure;

FIG. 7 is a schematic diagram of another structure and points randomlydistributed therein provided by the present disclosure;

FIG. 8 is a schematic diagram of a displacement field of the structurein FIG. 7 in x direction;

FIG. 9 is a schematic diagram of a stress field of the structure in FIG.7 in x direction;

FIG. 10 is a schematic diagram of a plate with a crack and points evenlydistributed therein provided by the present disclosure;

FIG. 11 is a schematic diagram of support domains of points P₁ and P₅respectively located at two sides of the crack and free boundariestherebetween provided by the present disclosure;

FIG. 12 is a schematic diagram of a displacement field of the plate witha crack in FIG. 10 in x direction; and

FIG. 13 is a schematic diagram of a stress field of the plate with acrack in FIG. 10 in y direction.

DETAILED DESCRIPTION

The present disclosure will be further described below with reference tothe accompanying drawings and specific embodiments. It should be notedthat, provided that there is no conflict, new embodiments can be formedby arbitrarily combining various embodiments or various technicalfeatures described below.

A meshless method for solid mechanics simulation specifically includesthe following three parts:

Part 1: Constructing a trial function and test a function.

For a solid structure of a given shape shown in FIG. 1, several pointsare distributed inside the structure and on its boundary. As shown inFIG. 2, with these points, the structure is divided into a set ofarbitrarily shaped subdomains, and each subdomain contains one point.The present disclosure does not limit the division of subdomains, merelyrequiring each subdomain to contain one point. For example, a VoronoiDiagram method can be adopted herein.

The trial function and the test function are constructed based on thesedistributed points and their subdomains, which specifically includes thefollowing steps.

Step A11: First, a polynomial trial function u^(h) (namely, adisplacement field, which refers to a displacement of the structure frombefore to after deformation) is defined for each subdomain. This trialfunction may be a first, second, third, or higher-order polynomialfunction.

For example, as shown in the following formula (1), nth-order Taylorexpansion is performed on the trial function at point P₀:

$\begin{matrix}{{{u^{h}\left( {x,y} \right)} = {\begin{bmatrix}u_{x}^{h} \\u_{y}^{h}\end{bmatrix} = \begin{bmatrix}{{u_{x}^{0} + \frac{\partial u_{x}}{\partial x}}❘_{P_{0}}{{\left( {x - x_{0}} \right) + \frac{\partial u_{x}}{\partial y}}❘_{P_{0}}{\left( {y - y_{0}} \right) + \ldots +}}} \\{\frac{1{\partial^{n}u_{x}}}{{n!}{\partial y^{n}}}\left( {y - y_{0}} \right)^{n}} \\{{u_{y}^{0} + \frac{\partial u_{y}}{\partial x}}❘_{P_{0}}{{\left( {x - x_{0}} \right) + \frac{\partial u_{y}}{\partial y}}❘_{P_{0}}{\left( {y - y_{0}} \right) + \ldots +}}} \\{\frac{1{\partial^{n}u_{y}}}{{n!}{\partial y^{n}}}\left( {y - y_{0}} \right)^{n}}\end{bmatrix}}},} & (1)\end{matrix}$

where [u_(x) ⁰ u_(y) ⁰]^(T) is value of the displacement u^(h) at thepoint P₀, and

$a = {\left\lbrack {\frac{\partial u_{x}}{\partial x}❘_{P_{0}}{\frac{\partial u_{x}}{\partial y}❘_{P_{0}}{{\ldots\frac{\partial^{n}u_{x}}{\partial y^{n}}\frac{\partial u_{y}}{\partial x}}❘_{P_{0}}{\frac{\partial u_{y}}{\partial y}❘_{P_{0}}{\ldots\frac{\partial^{n}u_{y}}{\partial y^{n}}}}}}} \right\rbrack ❘_{P_{0}}}$is a displacement gradient from the first order to the nth order at thepoint P₀.

When n=1, that is, when the trial function is a first-order function,the trial function u^(h) of the subdomain E₀ which contains the point P₀may be expressed as follows:

$\begin{matrix}{{{u^{h}\left( {x,y} \right)} = {\left\lbrack \text{⁠}\begin{matrix}u_{x}^{h} \\u_{y}^{h}\end{matrix} \right\rbrack = {{\begin{bmatrix}{{u_{x}^{0} + \frac{\partial u_{x}}{\partial x}}❘_{P_{0}}{{\left( {x - x_{0}} \right) + \frac{\partial u_{x}}{\partial y}}❘_{P_{0}}\left( {y - y_{0}} \right)}} \\{{u_{y}^{0} + \frac{\partial u_{y}}{\partial x}}❘_{P_{0}}{{\left( {x - x_{0}} \right) + \frac{\partial u_{y}}{\partial y}}❘_{P_{0}}\left( {y - y_{0}} \right)}}\end{bmatrix}\left( {x,y} \right)} \in E_{0}}}},} & (1)\end{matrix}$

where u^(h) (x, y) is value of u^(h) at point (x, y), (x₀, y₀) is thecoordinates of the point P₀, [u_(x) ⁰ u_(y) ⁰]^(T) is the value of u^(h)at the point P₀, and the displacement gradient at the point P₀ is

$a = \left\lbrack {\begin{matrix}\frac{\partial u_{x}}{\partial x} & \frac{\partial u_{x}}{\partial y} & \frac{\partial u_{y}}{\partial x} & \left. \frac{\partial u_{y}}{\partial y} \right\rbrack^{T}\end{matrix}❘_{P_{0}}} \right.$

In order to express the trial function with point displacements, it isrequired to express the displacement gradient with the pointdisplacements. The point displacements refer to displacements by whichthe point displaces, for example, displacements of black points in FIGS.2 and 3.

To express the displacement gradient with the point displacements, firstit is required to set the support domain for each point. In the presentdisclosure, the support domain may be defined in many manners, forexample, the support domain of each point is defined to contain allsubdomains having points adjacent to the corresponding point locatedtherein, as shown in FIG. 3. Definitely, the following definition mayalso be adopted: defining the support domain of each point to containall subdomains having points therein which are located in a circle of aparticular size and with the corresponding point at the center. Inaddition, if a crack exists between neighboring subdomains, the pointlocated on one side of the crack is not contained in the support domainof the point on the other side of the crack.

Step A12: The support domain of each subdomain is determined, and aweighted discrete L² norm J is defined, where

$\begin{matrix}{{J = {\left( {{Aa} + u_{0} - u_{m}} \right)^{T}{W\left( {{Aa} + u_{0} - u_{m}} \right)}}},} & (2)\end{matrix}$ $\begin{matrix}{A = \begin{bmatrix}{x_{1} - x_{0}} & {y_{1} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{1} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\{0} & 0 & {\ldots} & 0 & {x_{1} - x_{0}} & {y_{1} - y_{0}} & {\ldots} & {\frac{1}{n!}\left( {y_{1} - y_{0}} \right)^{n}} \\{x_{2} - x_{0}} & {y_{2} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{2} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & {x_{2} - x_{0}} & {y_{2} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{2} - y_{0}} \right)^{n}} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\{x_{m} - x_{0}} & {y_{m} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{m} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & {x_{m} - x_{0}} & {y_{m} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{m} - y_{0}} \right)^{n}}\end{bmatrix}} & \end{matrix}$ $\begin{matrix}{{u_{0} = \left( \begin{bmatrix}{u_{x}^{0}} & u_{y}^{0} & u_{x}^{0} & u_{y}^{0} & \ldots & u_{x}^{0} & u_{y}^{0}\end{bmatrix}_{1 \times 2m} \right)^{T}},} & \end{matrix}$ $\begin{matrix}{{u_{m} = \begin{bmatrix}{u_{x}^{1}} & u_{y}^{1} & u_{x}^{2} & u_{y}^{2} & \ldots & u_{x}^{m} & u_{y}^{m}\end{bmatrix}^{T}},} & \end{matrix}$ $\begin{matrix}{{W = \begin{bmatrix}w_{x}^{1} & 0 & 0 & \ldots & 0 \\0 & w_{y}^{1} & \ddots & \ddots & \vdots \\0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & w_{x}^{m} & 0 \\0 & \ldots & 0 & 0 & w_{y}^{m}\end{bmatrix}},} & \end{matrix}$

(x_(i), y_(i)) is coordinates of the point P_(i), [u_(x) ^(i) u_(y)^(i)]^(T) is value of the displacement u^(h) at the point P_(i), and[w_(x) ^(i) w_(y) ^(i)]^(T) is value of a weight function at the pointP_(i), (i=1, 2, 3, . . . m).

Step A13: The displacement gradient a is expressed with the pointdisplacements according to the formula (2) and stationary conditions ofthe norm J, that is:

$\begin{matrix}{{a = {Cu}_{E}},} & (3)\end{matrix}$ $\begin{matrix}{{{{wherein}C} = {\left( {A^{T}{WA}} \right)^{- 1}A^{T}{W\begin{bmatrix}I_{1} & I_{2}\end{bmatrix}}}},} & \end{matrix}$ $\begin{matrix}{{u_{E} = \begin{bmatrix}{u_{x}^{0}} & u_{y}^{0} & u_{x}^{1} & u_{y}^{1} & \ldots & u_{x}^{m} & u_{y}^{m}\end{bmatrix}^{T}},} & \end{matrix}$ $\begin{matrix}{{I_{1} = \begin{bmatrix}{- 1} & 0 \\0 & {- 1} \\{- 1} & 0 \\0 & {- 1} \\ \vdots & \vdots \\{- 1} & 0 \\0 & {- 1}\end{bmatrix}_{2m \times 2}},{{{{and}I_{2}} = \begin{bmatrix}1 & 0 & \ldots & 0 \\0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\0 & \ldots & 0 & 1\end{bmatrix}_{2m \times 2m}};}} & \end{matrix}$

Step A14: The following expression of the trial function in eachsubdomain is obtained according to formulae (3) and (1):

$\begin{matrix}{{u^{h} = {Nu}_{E}},} & (4)\end{matrix}$ $\begin{matrix}{N = {{\begin{bmatrix}{x - x_{0}} & {y - y_{0}} & 0 & 0 \\0 & 0 & {x - x_{0}} & {y - y_{0}}\end{bmatrix}C} + I_{3}}} & \end{matrix}$ $\begin{matrix}{where} & \end{matrix}$ $\begin{matrix}{{= \begin{bmatrix}N_{0} & 0 & N_{1} & 0 & \ldots & N_{m} & 0 \\0 & N_{0} & 0 & N_{1} & \ldots & 0 & N_{m}\end{bmatrix}_{2 \times {({{2m} + 2})}}},} & \end{matrix}$and N is referred to as a shape function herein; and

$I_{3} = {\begin{bmatrix}1 & 0 & 0 & \ldots & 0 \\0 & 1 & 0 & \ldots & 0\end{bmatrix}_{2 \times {({{2m} + 2})}}.}$

The test function and the trial function are stipulated to haveexpressions of an identical form, that is, they have the same shapefunction. Therefore, the following expression of the test function v canbe obtained by the same method: v=Nv_(E)(5).

Thus, according to the relationship between the displacement field andstrain and stress fields, the strain field ε^(h) and the stress fieldσ^(h) can be expressed with the trial functions and the test functionsrespectively as follows:

${\varepsilon^{h} = {\begin{bmatrix}\varepsilon_{x}^{h} \\\varepsilon_{y}^{h} \\\gamma_{xy}^{h}\end{bmatrix} = {{\begin{bmatrix}\frac{\partial}{\partial x} & 0 \\0 & \frac{\partial}{\partial y} \\\frac{\partial}{\partial y} & \frac{\partial}{\partial x}\end{bmatrix}u^{h}} = {{\begin{bmatrix}1 & 0 & 0 & 0 \\0 & 0 & 0 & 1 \\0 & 1 & 1 & 0\end{bmatrix}{Cu}_{E}} = {Bu}_{E}}}}}{{\sigma^{h} = {\begin{bmatrix}\sigma_{x}^{h} \\\sigma_{y}^{h} \\\sigma_{xy}^{h}\end{bmatrix} = {{D\begin{bmatrix}\varepsilon_{x} \\\varepsilon_{y} \\\gamma_{xy}\end{bmatrix}} = {DBu}_{E^{-}}}}},}$

where D is a material stiffness matrix and is related to a materialproperty of the structure, and the material property may be isotropy,orthotropy, anisotropy or another property.

Part 2: Introducing numerical flux to correct the Galerkin weak form.

Expressions of trial and test functions of the whole structure can beobtained according to the expressions of the trial and test functions ofeach subdomain.

The trial function and the test function of each subdomain arecontinuous within the subdomain. However, for the whole structure,because the subdomains have different displacement fields on a commonboundary, the trial function and the test function are discontinuous inthe whole structure.

For example, FIG. 5 is a schematic diagram of the trial function, and itcan be seen from the figure that the trial function in the wholestructure is discontinuous. FIG. 4 is a schematic diagram of a shapefunction of a point, and likewise the shape function is alsodiscontinuous.

The present disclosure is developed to be applied in the existing EFGmethod, and the existing Galerkin weak form has the following formula:

$\begin{matrix}{{\sum\limits_{E \in \Omega}{\int_{E}{{\sigma_{ij}(u)}{\varepsilon_{ij}(v)}d\Omega}}} = {{\sum\limits_{E \in \Omega}{\int_{E}{f_{i}v_{i}d\Omega}}} + {\sum\limits_{e \in \Gamma_{i}}{\int_{e}{{\overset{\_}{t}}_{l}v_{i}d{\Gamma.}}}}}} & (6)\end{matrix}$

However, because the trial and test functions of neighboring subdomainsgiven by the present disclosure are discontinuous, mathematicalinconsistency will be caused if the trial and test functions of eachsubdomain calculated in the present disclosure are directly used in theformula of the Galerkin weak form, finally resulting in instability ofcalculation results.

To resolve the instability of the calculation results that is caused bythe discontinuity of the trial and test functions, the presentdisclosure introduces numerical flux corrections to correct the existingGalerkin weak form.

The numerical flux corrections in the present disclosure may be ofmultiple types. Due to being capable of limiting discontinuity of afunction on an internal boundary, the numerical flux correction can beintroduced herein to resolve the instability of numerical results causedby the discontinuity of the trial function.

This embodiment uses interior penalty numerical flux corrections as anexample to correct the existing Galerkin weak form. The followingformula (7) indicates the Galerkin weak form with the interior penalty(IP) numerical flux corrections. In the formula (7), traction boundaryconditions and displacement continuity conditions of the internalboundary are imposed weakly, while displacement boundary conditions arerequired to be imposed strongly:

$\begin{matrix}{{{\sum\limits_{E \in \Omega}{\int_{E}{{\sigma_{ij}(u)}{\varepsilon_{ij}(v)}d\Omega}}} - {\sum\limits_{e \in \Gamma_{k}}{\int_{e}{{\left\{ {{\sigma_{ij}(u)}n_{j}^{e}} \right\}\left\lbrack v_{i} \right\rbrack}d\Gamma}}} - {\sum\limits_{e \in \Gamma_{k}}{\int_{e}{{\left\{ {{\sigma_{ij}(v)}n_{j}^{e}} \right\}\left\lbrack u_{i} \right\rbrack}d\Gamma}}} + {\sum\limits_{e \in \Gamma_{k}}{\frac{\eta}{h_{\in}}{\int_{e}{{\left\lbrack u_{i} \right\rbrack\left\lbrack v_{i} \right\rbrack}d\Gamma}}}}} = {{\sum\limits_{E \in \Omega}{\int_{E}{f_{i}v_{i}d\Omega}}} + {\sum\limits_{e \in \Gamma_{k}}{\int_{e}{{\overset{\_}{t}}_{l}v_{i}d{\Gamma.}}}}}} & (7)\end{matrix}$

In the formula, Γ_(u) is the displacement boundary, Γ_(t) is thetraction boundary, { } is an average operator, and [ ] is a jumpoperator. When e is the internal boundary (a common boundary betweensubdomains E₁ and E₂), n_(j) ^(e) indicates a unit normal vectorpointing from

${E_{1}{to}E_{2}},{\lbrack w\rbrack = {w❘_{e}^{E_{1}}{{- w}❘_{e}^{E_{2}}}}},{{{and}\left\{ w \right\}} = {\frac{1}{2}{\left( {w❘_{e}^{E_{1}}{{+ w}❘_{e}^{E_{2}}}} \right).}}}$

It can be known from the above that, by comparison between the formula(7) indicating the Galerkin weak form with the IP numerical fluxcorrections provided by the present disclosure and the formula (6)indicating the traditional Galerkin weak form, their difference lies inthat the former has a jump and an average operator on the commonboundary, thus resolving the inconsistency of numerical results causedby the discontinuity of the trial and test functions.

Part 3: Substituting the expressions of the trial and test functions ineach subdomain into the Galerkin weak form with the numerical fluxcorrections, to construct a global stiffness matrix of the structure anda system of algebraic equations; calculating the displacement of eachpoint in the structure according to the system of algebraic equations;and afterwards, calculating deformation and stress of each position inthe structure according to the displacement of each point in thestructure, completing simulation calculation for the structure.

Step B1: First, the trial and test functions in each subdomain aresubstituted into the Galerkin weak form with the numerical fluxcorrections, that is, the trial and test functions in each subdomain aresubstituted into the formula (7), to obtain a point stiffness matrix ofeach point and an internal-boundary stiffness matrix.

Where the point stiffness matrix is K_(E)=∫_(E)B^(T)DBdΩ, theinternal-boundary stiffness matrix is

${K_{h} = {{- {\int_{e}{\left\lbrack N^{T} \right\rbrack\left\{ {n_{e}{DB}} \right\} d\Gamma}}} - {\int_{e}{{\left\{ {B^{T}{Dn}_{e}^{T}} \right\}\lbrack N\rbrack}d\Gamma}} + {\frac{\eta}{h_{e}}{\int_{e}{{\left\lbrack N^{T} \right\rbrack\lbrack N\rbrack}d\Gamma}}}}},$and e is a common boundary between subdomains.

Step B2: The global stiffness matrix K of the structure is obtained byassembling the point stiffness matrix of each point and theinternal-boundary stiffness matrix, to obtain a system of algebraicequations in a matrix form:Kq=Q  (8),where q is a point displacement vector to be solved, Q is a load vector,and K is the global stiffness matrix obtained by assembling all thepoint stiffness matrices and the internal-boundary stiffness matrices.

Step B3: The system of algebraic equations shown by the formula (8) issolved according to the global stiffness matrix of the structure and acurrent load vector, to obtain a displacement vector of each point inthe structure, namely, a displacement of each point.

Step B4: Displacements of each position in the structure are obtainedaccording to the displacement of each point, and strain and stress ofeach position in the structure is obtained according to the relationshipbetween the displacement and the stress.

By the present disclosure, analysis results regarding structuraldeformation and stress can be obtained. When the analysis results areapplied in a structural design, it can be determined whether designs ofobject structures such as mechanical and civil engineering structuresare safe and reasonable.

The load, which may be a single load condition or a load history, refersto a force imposed on a structure in a working status. The samestructure has different deformations and stresses under the effect ofdifferent loads. If the present disclosure is applied in a structuresubjected to a load history in a working status, this load history maybe divided into m load steps. In each load step, deformation anddisplacement are calculated and crack generation and propagation aresimulated.

The present disclosure further provides a calculation process ofstructural strain and stress with reference to specific drawings.

Embodiment 1

For a structure shown in FIG. 7, the following steps are performed.

Step C11: The structure is divided into several subdomains according topoints randomly distributed in the structure.

Step C12: the support domain of each point is determined, andexpressions of trial and test functions in each subdomain in thestructure are calculated by using a weighted least square method.

The expression of the trial function is u^(h)=Nu_(E), N is a shapefunction, and u_(E) is a point displacement.

The expression of the test function is v=Nv_(E), which has the sameshape function N as the trial function.

Step C13: numerical flux corrections are introduced to obtain theGalerkin weak form with numerical flux corrections.

Step C14: The trial and test functions in each subdomain are substitutedinto the Galerkin weak form with the numerical flux corrections, namely,the formula (7); and then numerical integration is performed by means ofGaussian integration to obtain a point stiffness matrix K_(E) and aninternal-boundary stiffness matrix K_(h). Because the trial and testfunctions are both polynomial functions, by use of Gaussian integrationto carry out numerical integration, the point stiffness matrix K_(E) andthe internal-boundary stiffness matrix K_(h) can be rapidly andaccurately calculated.

Step C15: A global stiffness matrix K of the structure is obtained byassembling all of the point stiffness matrices K_(E) and theinternal-boundary stiffness matrices K_(h), and the following system ofalgebraic equations is obtained according to a given load vector: Kq=Q,where Q is the given load vector and q is a point displacement vector.

The point displacement vector q is calculated by strongly imposingdisplacement boundary conditions, and then is substituted back into theexpression of the trial function, to obtain a displacement field of thestructure.

A stress field of the structure is obtained according to therelationship between the displacement field and the stress field.

For example, a displacement field of the structure in FIG. 7 in xdirection is shown in FIG. 8, and a stress field thereof in x directionis shown in FIG. 9. By comparison between the foregoing simulationresults and exact solution in theory, the relative error of thedisplacement field is 0.051% and that of the stress field is 2.064%.

Embodiment 2

For a plate with a crack in FIG. 10, the following steps are performed:

Step C21: The structure is divided into several subdomains by means ofthe identical subdomain division method.

Step C22: the support domain of each point is determined. Because thereis a crack in the structure, it is assumed that points located on oneside of the crack are not contained in the support domain of points onthe other side of the crack. As shown in FIG. 11, the support domain ofthe point P₁ does not contain the point P₅, and the support domain ofthe point P₅ does not contain the point P₁.

Moreover, boundaries that coincide with the crack are set astraction-free boundaries. As shown in FIG. 11, boundaries Γ₂₆, Γ₁₅, andΓ₄₇ are all traction-free boundaries.

Step C23: The following expression of the trial function in eachsubdomain is calculated by using a weighted least square method:u^(h)=Nu_(E), where N is a shape function and u_(E) is a pointdisplacement.

An expression of the test function is v=Nv_(E), which has the same shapefunction as the trial function.

Step C24: The trial and test functions in each subdomain are substitutedinto the Galerkin weak form with the numerical flux corrections, namely,the formula (7); and then numerical integration is performed by means ofGaussian integration to obtain a point stiffness matrix K_(E) and aninternal-boundary stiffness matrix K_(h).

A global stiffness matrix K is obtained by assembling all of the pointstiffness matrices K_(E) and the internal-boundary stiffness matricesK_(h), and the following system of algebraic equations is obtainedaccording to a current load vector: Kq=Q, where Q is the current loadvector and q is a point displacement vector.

The point displacement vector q is calculated by strongly imposingdisplacement boundary conditions, and then is substituted back into theexpression of the trial function, to obtain the displacement field ofthe structure.

The stress field of the structure is obtained according to therelationship between the displacement field and the stress field.

A displacement field of the plate with a crack in FIG. 10 in x directionis shown in FIG. 12, and a stress field thereof in y direction is shownin FIG. 13. An error between a stress intensity factor calculated byusing a J integral in the present disclosure and an exact value intheory is 0.83%.

In addition, the structure in a working status may not only be subjectedto a static load, but may also experience a load history. The loadhistory herein includes a series of loads in chronological order. Inthis case, the load history may be divided into m load steps, and themethod is performed in each load step. Afterwards, according to thecalculated deformation and stress of each position in the structure, itis determined, according to a certain criterion, whether there is acrack between neighboring subdomains in the structure. Thus, rupture andfragmentation of the structure under the effect of a certain loadhistory can be simulated.

For the static load which refers to that the structure is subjected toonly one load in a working status, by using the meshless method forsolid mechanics simulation provided by the present disclosure,deformation and stress of each position in the structure can beobtained, completing simulation for the structure.

For the load history which refers to that the structure is subjected toa series of loads in chronological order in a working status, the stepsin the meshless method for solid mechanics simulation provided by thepresent disclosure are performed under each of the chronologicallyarranged loads successively, to obtain deformation and stress of eachposition in the structure; and then it is determined, according to acertain criterion, whether there is a crack between neighboringsubdomains in the structure under a corresponding load. Thus, ruptureand fragmentation of the structure under the effect of the load historycan be simulated, as shown in FIG. 6.

The present disclosure achieves the following technical effects comparedto the prior art:

1. The present disclosure provides polynomial and discontinuous trialand test functions.

The present disclosure divides a structure into several subdomains bymeans of novel point-based interpolation, and then calculates the trialfunction and the test function of each subdomain. Thus, the resultanttrial function and test function are both polynomial functions in eachsubdomain. Because the trial function and the test function havedifferent values on the common boundary between neighboring subdomains,the trial and test functions of each subdomain are discontinuous amongdifferent subdomains.

2. The present disclosure further corrects the existing Galerkin weakform. Inconsistency of numerical results that is caused by discontinuityof the trial and test functions provided by the present disclosure isresolved by introducing numerical flux corrections.

3. Because the trial and test functions are polynomial functions in thepresent disclosure, a stiffness matrix of the structure is calculated bymeans of Gaussian integration or an analytical method, so that integralscan be computed rapidly, easily, and accurately.

4. The present disclosure makes for a symmetric, sparse and positivedefinite global stiffness matrix, thus being applicable to large-scalestructural simulation.

5. Due to the discontinuity of the trial and test functions, the presentdisclosure can easily cut off the connection between points, and furthercan easily and efficiently simulate crack generation and propagation,and rupture and fragmentation of a structure according to a certaincriterion for cracking.

The present disclosure further provides an electronic device, whichincludes a memory, a processor, and a computer program stored in thememory and executable in the processor. The processor implements thesteps in the meshless method for solid mechanics simulation described inthe specification when executing the program.

The present disclosure further provides a computer-readable storagemedium which stores a computer program. When executed by a processor,the computer program implements the steps in the meshless method forsolid mechanics simulation described in the specification.

The above merely describes preferred embodiments of the presentdisclosure, and is not intended to limit the scope of protection of thepresent disclosure. Any non-substantial changes and replacements made onthe basis of the present disclosure by persons skilled in the art allfall within the scope of protection of the present disclosure.

What is claimed is:
 1. A meshless method for solid mechanics simulation,configured for calculating deformation and stress of a solid structureof a given shape, comprising the following steps: a step of constructinga trial function and a test function executed by a processor of anelectronic device, wherein obtaining expressions of the trial functionand the test function with point displacements in each subdomain of thesolid structure of a given shape according to points distributed in thesolid structure of a given shape, wherein the trial function and thetest function of neighboring subdomains of the solid structure of agiven shape are discontinuous; a step of developing a Galerkin weak formexecuted by the processor of the electronic device, wherein introducingnumerical flux corrections to obtain the Galerkin weak form with thenumerical flux corrections; a step of generating and solving a system ofalgebraic equations executed by the processor of the electronic device,wherein constructing a global stiffness matrix of the solid structure ofa given shape and the system of algebraic equations according to theexpressions of the trial function and test function in each subdomain ofthe solid structure of a given shape and the Galerkin weak form with thenumerical flux corrections, and solving the system of algebraicequations to obtain a displacement of each point; and a step ofcalculating structural deformation and stress executed by the processorof the electronic device, wherein calculating the deformation and thestress of each position in the solid structure of a given shapeaccording to the displacement of each point, to complete simulation andanalysis of the solid structure of a given shape; wherein, the methodfurther comprising: a simulating step, executed by the processor of theelectronic device, for according to a sequence of loads in chronologicalorder, performing the step of constructing a trial function and a testfunction, the step of developing a Galerkin weak form, the step ofgenerating and solving a system of algebraic equations, and the step ofcalculating structural deformation and stress under each loadsuccessively, to obtain the deformation and stress of each position inthe solid structure of a given shape; and determining, according to acertain criterion, whether there is a crack between neighboringsubdomains in the solid structure of a given shape, to simulate ruptureand fragmentation of the solid structure of a given shape under theeffect of the sequence of loads in chronological order.
 2. The meshlessmethod for solid mechanics simulation according to claim 1, wherein eachsubdomain contains one point.
 3. The meshless method for solid mechanicssimulation according to claim 2, wherein a process of constructing thetrial function and the test function is as follows: step S11: firstdefining the trial function in each subdomain and recording the trialfunction as u^(h), wherein the trial function is an nth-order polynomialfunction, and n is a natural number greater than or equal to 1 and oftenequal to 1; and performing nth-order Taylor expansion on the trialfunction at a point P₀, as shown in the following formula (1):$\begin{matrix}{{{u^{h}\left( {x,y} \right)} = {\begin{bmatrix}u_{x}^{h} \\u_{y}^{h}\end{bmatrix} = \begin{bmatrix}{{u_{x}^{0} + \frac{\partial u_{x}}{\partial x}}❘_{P_{0}}{{\left( {x - x_{0}} \right) + \frac{\partial u_{x}}{\partial y}}❘_{P_{0}}{\left( {y - y_{0}} \right) + \ldots +}}} \\{\frac{1{\partial^{n}u_{x}}}{{n!}{\partial y^{n}}}\left( {y - y_{0}} \right)^{n}} \\{{u_{y}^{0} + \frac{\partial u_{y}}{\partial x}}❘_{P_{0}}{{\left( {x - x_{0}} \right) + \frac{\partial u_{y}}{\partial y}}❘_{P_{0}}{\left( {y - y_{0}} \right) + \ldots +}}} \\{\frac{1{\partial^{n}u_{y}}}{{n!}{\partial y^{n}}}\left( {y - y_{0}} \right)^{n}}\end{bmatrix}}},} & (1)\end{matrix}$ wherein [u_(x) ⁰ u_(y) ⁰]^(T) is value of the displacementu^(h) at the point P₀, and$a = {\left\lbrack {\frac{\partial u_{x}}{\partial x}❘_{P_{0}}{\frac{\partial u_{x}}{\partial y}❘_{P_{0}}{{\ldots\frac{\partial^{n}u_{x}}{\partial y^{n}}\frac{\partial u_{y}}{\partial x}}❘_{P_{0}}{\frac{\partial u_{y}}{\partial x}❘_{P_{0}}{\ldots\frac{\partial^{n}u_{y}}{\partial y^{n}}}}}}} \right\rbrack ❘_{P_{0}}}$is a displacement gradient from the first order to the nth order at thepoint P₀, where n is the order of Taylor expansion; step S12:determining support domain of each subdomain, and defining a weighteddiscrete L² norm J, wherein $\begin{matrix}{{J = {\left( {{Aa} + u_{0} - u_{m}} \right)^{T}{W\left( {{Aa} + u_{0} - u_{m}} \right)}}},} & (2)\end{matrix}$ $\begin{matrix}{A = \begin{bmatrix}{x_{1} - x_{0}} & {y_{1} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{1} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\{0} & 0 & {\ldots} & 0 & {x_{1} - x_{0}} & {y_{1} - y_{0}} & {\ldots} & {\frac{1}{n!}\left( {y_{1} - y_{0}} \right)^{n}} \\{x_{2} - x_{0}} & {y_{2} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{2} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & {x_{2} - x_{0}} & {y_{2} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{2} - y_{0}} \right)^{n}} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\{x_{m} - x_{0}} & {y_{m} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{m} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & {x_{m} - x_{0}} & {y_{m} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{m} - y_{0}} \right)^{n}}\end{bmatrix}} & \end{matrix}$ $\begin{matrix}{{u_{0} = \left( \begin{bmatrix}{u_{x}^{0}} & u_{y}^{0} & u_{x}^{0} & u_{y}^{0} & \ldots & u_{x}^{0} & u_{y}^{0}\end{bmatrix}_{1 \times 2m} \right)^{T}},} & \end{matrix}$ $\begin{matrix}{{u_{m} = \begin{bmatrix}{u_{x}^{1}} & u_{y}^{1} & u_{x}^{2} & u_{y}^{2} & \ldots & u_{x}^{m} & u_{y}^{m}\end{bmatrix}^{T}},} & \end{matrix}$ $\begin{matrix}{{W = \begin{bmatrix}w_{x}^{1} & 0 & 0 & \ldots & 0 \\0 & w_{y}^{1} & \ddots & \ddots & \vdots \\0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & w_{x}^{m} & 0 \\0 & \ldots & 0 & 0 & w_{y}^{m}\end{bmatrix}},} & \end{matrix}$ (x_(i), y_(i)) is coordinates of a point P_(i), [u_(x)^(i) u_(y) ^(i)]^(T) is value of the displacement u^(h) at the pointP_(i), and [w_(x) ^(i) w_(y) ^(i)]^(T); is value of a weight function atthe point P_(i), (i=1, 2, 3, . . . m); step S13: expressing thedisplacement gradient a with the point displacement according to theformula (2) and stationary conditions of the norm J, that is:$\begin{matrix}{{a = {Cu}_{E}},} & (3)\end{matrix}$ $\begin{matrix}{{{{wherein}C} = {\left( {A^{T}{WA}} \right)^{- 1}A^{T}{W\begin{bmatrix}I_{1} & I_{2}\end{bmatrix}}}},} & \end{matrix}$ $\begin{matrix}{{u_{E} = \begin{bmatrix}u_{x}^{0} & u_{y}^{0} & u_{x}^{1} & u_{y}^{1} & \ldots & u_{x}^{m} & u_{y}^{m}\end{bmatrix}^{T}},} & \end{matrix}$ $\begin{matrix}{{I_{1} = \begin{bmatrix}{- 1} & 0 \\0 & {- 1} \\{- 1} & 0 \\0 & {- 1} \\ \vdots & \vdots \\{- 1} & 0 \\0 & {- 1}\end{bmatrix}_{2m \times 2}},{{{{and}I_{2}} = \begin{bmatrix}1 & 0 & \ldots & 0 \\0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\0 & \ldots & 0 & 1\end{bmatrix}_{2m \times 2m}};}} & \end{matrix}$ step S14: obtaining the expression of the trial functionin each subdomain according to formulae (3) and (1); and u^(h) = Nu_(E)${ɛ^{h} = {\begin{bmatrix}ɛ_{11}^{h} \\ɛ_{22}^{h} \\{2ɛ_{12}^{h}}\end{bmatrix} = {{\begin{bmatrix}\frac{\partial}{\partial x_{1}} & 0 \\0 & \frac{\partial}{\partial x_{2}} \\\frac{\partial}{\partial x_{2}} & \frac{\partial}{\partial x_{1}}\end{bmatrix}u^{h}} = {Bu_{E}}}}},{\sigma^{h} = {\begin{bmatrix}\sigma_{11}^{h} \\\sigma_{22}^{h} \\\sigma_{12}^{h}\end{bmatrix} = {{Dɛ^{h}} = {DBu}_{E}}}}$ N is the shape function; D isthe stress-strain matrix; B is the derivative of shape function; ε^(h)is the strain and σ^(h) is the stress; step S15: obtaining theexpression of the test function by prescribing the test function in asame form as that of the trial function.
 4. The meshless method forsolid mechanics simulation according to claim 3, wherein a definition ofthe support domain comprises any one of the following manners: definingthe support domain of each point to contain all subdomains having pointsadjacent to the corresponding point located therein, or defining thesupport domain of each point to contain all subdomains having pointstherein which are located in a circle of a particular size and with thecorresponding point at a center; and if a crack exists between twoneighboring subdomains, the point located on one side of the crack isnot contained in the support domain of the point on the other side ofthe crack.
 5. The meshless method for solid mechanics simulationaccording to claim 1, wherein the numerical flux corrections compriseinterior penalty numerical flux corrections.
 6. The meshless method forsolid mechanics simulation according to claim 3, wherein the generatingand solving step of the system of algebraic equations further comprises:step S21: substituting the trial function and test function in eachsubdomain into the Galerkin weak form corrected by numerical fluxes, toobtain a point stiffness matrix and an internal-boundary stiffnessmatrix, wherein the point stiffness matrix is K_(E)=∫_(E)B^(T)DBdΩ, theinternal-boundary stiffness matrix is${K_{h} = {{- {\int_{e}{\left\lbrack N^{T} \right\rbrack\left\{ {n_{e}{DB}} \right\} d\Gamma}}} - {\int_{e}{{\left\{ {B^{T}{Dn}_{e}^{T}} \right\}\lbrack N\rbrack}d\Gamma}} + {\frac{\eta}{h_{e}}{\int_{e}{{\left\lbrack N^{T} \right\rbrack\lbrack N\rbrack}d\Gamma_{t}}}}}},$and e is a common boundary between neighboring subdomains; step S22:obtaining a global stiffness matrix K by assembling all the pointstiffness matrices and the internal-boundary stiffness matrices, toobtain the following system of algebraic equations in a matrix form:Kq=Q, wherein q is a point displacement vector, and Q is a load vector;and step S23: solving the system of algebraic equations according to theglobal stiffness matrix and the load vector to obtain the pointdisplacement vector of the solid structure of a given shape, thusfinally obtaining the displacement of each position in the solidstructure of a given shape.
 7. An electronic device, comprising amemory, a processor, and a computer program stored in the memory andexecutable in the processor, wherein the processor implements thefollowing steps when executing the program a step of constructing atrial function and a test function, wherein obtaining expressions of thetrial function and the test function with point displacements in eachsubdomain according to points distributed in a solid structure of agiven shape, wherein the trial function and the test function ofneighboring subdomains of the solid structure of a given shape arediscontinuous; a step of developing a Galerkin weak form, whereinintroducing numerical flux corrections to obtain the Galerkin weak formwith the numerical flux corrections; a step of generating and solving asystem of algebraic equations, wherein constructing a global stiffnessmatrix of the solid structure of a given shape and the system ofalgebraic equations according to the expressions of the trial functionand test function in each subdomain and the Galerkin weak form with thenumerical flux corrections, and solving the system of algebraicequations to obtain a displacement of each point; and a step ofcalculating structural deformation and stress, wherein calculating thedeformation and the stress of each position in the solid structure of agiven shape according to the displacement of each point, to completesimulation and analysis of the structure; a simulating step foraccording to a sequence of loads in chronological order, performing thestep of constructing a trial function and a test function, the step ofdeveloping a Galerkin weak form, the step of generating and solving asystem of algebraic equations, and the step of calculating structuraldeformation and stress under each load successively, to obtain thedeformation and stress of each position in the solid structure of agiven shape; and determining, according to a certain criterion, whetherthere is a crack between neighboring subdomains in the solid structureof a given shape, to simulate rupture and fragmentation of the solidstructure of a given shape under the effect of the sequence of loads inchronological order.
 8. The electronic device according to claim 7,wherein each subdomain contains one point.
 9. The electronic deviceaccording to claim 8, wherein a process of constructing the trialfunction and the test function is as follows: step S11: first definingthe trial function in each subdomain and recording the trial function asu^(h), wherein the trial function is an nth-order polynomial function,and n is a natural number greater than or equal to 1 and often equal to1; and performing nth-order Taylor expansion of the trial function at apoint P₀, as shown in the following formula (1): $\begin{matrix}{{{u^{h}\left( {x,y} \right)} = {\begin{bmatrix}u_{x}^{h} \\u_{y}^{h}\end{bmatrix} = \begin{bmatrix}{{u_{x}^{0} + \frac{\partial u_{x}}{\partial x}}❘_{P_{0}}{{\left( {x - x_{0}} \right) + \frac{\partial u_{x}}{\partial y}}❘_{P_{0}}{\left( {y - y_{0}} \right) + \ldots +}}} \\{\frac{1{\partial^{n}u_{x}}}{{n!}{\partial y^{n}}}\left( {y - y_{0}} \right)^{n}} \\{{u_{y}^{0} + \frac{\partial u_{y}}{\partial x}}❘_{P_{0}}{{\left( {x - x_{0}} \right) + \frac{\partial u_{y}}{\partial y}}❘_{P_{0}}{\left( {y - y_{0}} \right) + \ldots +}}} \\{\frac{1{\partial^{n}u_{y}}}{{n!}{\partial y^{n}}}\left( {y - y_{0}} \right)^{n}}\end{bmatrix}}},} & (1)\end{matrix}$ wherein [u_(x) ⁰ u_(y) ⁰]^(T) is value of the displacementu^(h) at the point P₀, and$a = {\left\lbrack {\frac{\partial u_{x}}{\partial x}❘_{P_{0}}{\frac{\partial u_{x}}{\partial y}❘_{P_{0}}{{\ldots\frac{\partial^{n}u_{x}}{\partial y^{n}}\frac{\partial u_{y}}{\partial x}}❘_{P_{0}}{\frac{\partial u_{y}}{\partial x}❘_{P_{0}}{\ldots\frac{\partial^{n}u_{y}}{\partial y^{n}}}}}}} \right\rbrack ❘_{P_{0}}}$is a displacement gradient from the first order to the nth order at thepoint P₀, where n is the order of Taylor expansion; step S12:determining support domain of each subdomain, and defining a weighteddiscrete L² norm J, wherein $\begin{matrix}{{J = {\left( {{Aa} + u_{0} - u_{m}} \right)^{T}{W\left( {{Aa} + u_{0} - u_{m}} \right)}}},} & (2)\end{matrix}$ $\begin{matrix}{A = \begin{bmatrix}{x_{1} - x_{0}} & {y_{1} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{1} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\{0} & 0 & {\ldots} & 0 & {x_{1} - x_{0}} & {y_{1} - y_{0}} & {\ldots} & {\frac{1}{n!}\left( {y_{1} - y_{0}} \right)^{n}} \\{x_{2} - x_{0}} & {y_{2} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{2} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & {x_{2} - x_{0}} & {y_{2} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{2} - y_{0}} \right)^{n}} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\{x_{m} - x_{0}} & {y_{m} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{m} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & {x_{m} - x_{0}} & {y_{m} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{m} - y_{0}} \right)^{n}}\end{bmatrix}} & \end{matrix}$ $\begin{matrix}{{u_{0} = \left( \begin{bmatrix}{u_{x}^{0}} & u_{y}^{0} & u_{x}^{0} & u_{y}^{0} & \ldots & u_{x}^{0} & u_{y}^{0}\end{bmatrix}_{1 \times 2m} \right)^{T}},} & \end{matrix}$ $\begin{matrix}{{u_{m} = \begin{bmatrix}{u_{x}^{1}} & u_{y}^{1} & u_{x}^{2} & u_{y}^{2} & \ldots & u_{x}^{m} & u_{y}^{m}\end{bmatrix}^{T}},} & \end{matrix}$ $\begin{matrix}{{W = \begin{bmatrix}w_{x}^{1} & 0 & 0 & \ldots & 0 \\0 & w_{y}^{1} & \ddots & \ddots & \vdots \\0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & w_{x}^{m} & 0 \\0 & \ldots & 0 & 0 & w_{y}^{m}\end{bmatrix}},} & \end{matrix}$ (x_(i), y_(i)) is coordinates of a point P_(i), [u_(x)^(i) u_(y) ^(i)]^(T) is value of the displacement u^(h) at the pointP_(i), and [w_(x) ^(i) w_(y) ^(i)]^(T) is value of a weight function atthe point P_(i), (i=1, 2, 3, . . . m); step S13: expressing thedisplacement gradient a with point displacements according to theformula (2) and stationary conditions of the norm J, that is:$\begin{matrix}{{a = {Cu}_{E}},} & (3)\end{matrix}$ $\begin{matrix}{{{{wherein}C} = {\left( {A^{T}{WA}} \right)^{- 1}A^{T}{W\begin{bmatrix}I_{1} & I_{2}\end{bmatrix}}}},} & \end{matrix}$ $\begin{matrix}{{u_{E} = \begin{bmatrix}u_{x}^{0} & u_{y}^{0} & u_{x}^{1} & u_{y}^{1} & \ldots & u_{x}^{m} & u_{y}^{m}\end{bmatrix}^{T}},} & \end{matrix}$ $\begin{matrix}{{I_{1} = \begin{bmatrix}{- 1} & 0 \\0 & {- 1} \\{- 1} & 0 \\0 & {- 1} \\ \vdots & \vdots \\{- 1} & 0 \\0 & {- 1}\end{bmatrix}_{2m \times 2}},{{{{and}I_{2}} = \begin{bmatrix}1 & 0 & \ldots & 0 \\0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\0 & \ldots & 0 & 1\end{bmatrix}_{2m \times 2m}};}} & \end{matrix}$ step S14: obtaining the expression of the trial functionin each subdomain according to formulae (3) and (1); and u^(h) = Nu_(E)${ɛ^{h} = {\begin{bmatrix}ɛ_{11}^{h} \\ɛ_{22}^{h} \\{2ɛ_{12}^{h}}\end{bmatrix} = {{\begin{bmatrix}\frac{\partial}{\partial x_{1}} & 0 \\0 & \frac{\partial}{\partial x_{2}} \\\frac{\partial}{\partial x_{2}} & \frac{\partial}{\partial x_{1}}\end{bmatrix}u^{h}} = {Bu_{E}}}}},{\sigma^{h} = {\begin{bmatrix}\sigma_{11}^{h} \\\sigma_{22}^{h} \\\sigma_{12}^{h}\end{bmatrix} = {{Dɛ^{h}} = {DBu}_{E}}}}$ N is the shape function; D isthe stress-strain matrix; B is the derivative of shape function; ε^(h)is the strain and σ^(h) is the stress; step S15: obtaining theexpression of the test function by prescribing the test function in asame form as that of the trial function.
 10. The electronic deviceaccording to claim 9, wherein a definition of the support domaincomprises any one of the following manners: defining the support domainof each point to contain all subdomains having points adjacent to thecorresponding point located therein, or defining the support domain ofeach point to contain all subdomains having points therein which arelocated in a circle of a particular size and with the correspondingpoint at a center; and if a crack exists between two neighboringsubdomains, the point located on one side of the crack is not containedin the support domain of the point on the other side of the crack. 11.The electronic device according to claim 9, wherein the step ofgenerating and solving the system of algebraic equations furthercomprises: step S21: substituting the trial function and test functionin each subdomain into the Galerkin weak form corrected by numericalfluxes, to obtain a point stiffness matrix and an internal-boundarystiffness matrix, wherein the point stiffness matrix isK_(E)=∫_(E)B^(T)DBdΩ, the internal-boundary stiffness matrix is${K_{h} = {{- {\int_{e}{\left\lbrack N^{T} \right\rbrack\left\{ {n_{e}{DB}} \right\} d\Gamma}}} - {\int_{e}{{\left\{ {B^{T}{Dn}_{e}^{T}} \right\}\lbrack N\rbrack}d\Gamma}} + {\frac{\eta}{h_{e}}{\int_{e}{{\left\lbrack N^{T} \right\rbrack\lbrack N\rbrack}d\Gamma_{t}}}}}},$and e is a common boundary between subdomains; step S22: obtaining aglobal stiffness matrix K by assembling all the point stiffness matricesand the internal-boundary stiffness matrices, to obtain the followingsystem of algebraic equations in a matrix form: Kq=Q, wherein q is apoint displacement vector, and Q is a load vector; and step S23: solvingthe system of algebraic equations according to the global stiffnessmatrix and the load vector to obtain the point displacement vector ofthe solid structure of a given shape, thus finally obtaining thedisplacement of each position in the solid structure of a given shape.12. The electronic device according to claim 7, wherein the numericalflux corrections comprise interior penalty numerical flux corrections.13. A non-transitory computer-readable storage medium which stores acomputer program, wherein when executed by a processor, the computerprogram implements the following steps a step of constructing a trialfunction and a test function, wherein obtaining expressions of the trialfunction and the test function with point displacements in eachsubdomain according to points distributed in a solid structure of agiven shape, wherein the trial function and the test function ofneighboring subdomains of the solid structure of a given shape arediscontinuous; a step of developing a Galerkin weak form, whereinintroducing numerical flux corrections to obtain the Galerkin weak formwith the numerical flux corrections; a step of generating and solving asystem of algebraic equations, wherein constructing a global stiffnessmatrix of the solid structure of a given shape and the system ofalgebraic equations according to the expressions of the trial functionand test function in each subdomain and the Galerkin weak form with thenumerical flux corrections, and solving the system of algebraicequations to obtain a displacement of each point; and a step ofcalculating structural deformation and stress, wherein calculating thedeformation and the stress of each position in the solid structure of agiven shape according to the displacement of each point, to completesimulation and analysis of the solid structure of a given shape;wherein, a simulating step for according to a sequence of loads inchronological order, performing the step of constructing a trialfunction and a test function, the step of developing a Galerkin weakform, the step of generating and solving a system of algebraicequations, and the step of calculating structural deformation and stressunder each load successively, to obtain the deformation and stress ofeach position in the solid structure of a given shape; and determining,according to a certain criterion, whether there is a crack betweenneighboring subdomains in the solid structure of a given shape, tosimulate rupture and fragmentation of the solid structure of a givenshape under the effect of the sequence of loads in chronological order.14. The non-transitory computer-readable storage medium according toclaim 13, wherein each subdomain contains one point.
 15. Thenon-transitory computer-readable storage medium according to claim 14,wherein a process of constructing the trial function and the testfunction is as follows: step S11: first defining the trial function ineach subdomain and recording the trial function as u^(h), wherein thetrial function is an nth-order polynomial function, and n is a naturalnumber greater than or equal to 1 and often equal to 1; and performingnth-order Taylor expansion of the trial function at a point P₀, as shownin the following formula (1): $\begin{matrix}{{{u^{h}\left( {x,y} \right)} = {\begin{bmatrix}u_{x}^{h} \\u_{y}^{h}\end{bmatrix} = \begin{bmatrix}{{u_{x}^{0} + \frac{\partial u_{x}}{\partial x}}❘_{P_{0}}{{\left( {x - x_{0}} \right) + \frac{\partial u_{x}}{\partial y}}❘_{P_{0}}{\left( {y - y_{0}} \right) + \ldots +}}} \\{\frac{1{\partial^{n}u_{x}}}{{n!}{\partial y^{n}}}\left( {y - y_{0}} \right)^{n}} \\{{u_{y}^{0} + \frac{\partial u_{y}}{\partial x}}❘_{P_{0}}{{\left( {x - x_{0}} \right) + \frac{\partial u_{y}}{\partial y}}❘_{P_{0}}{\left( {y - y_{0}} \right) + \ldots +}}} \\{\frac{1{\partial^{n}u_{y}}}{{n!}{\partial y^{n}}}\left( {y - y_{0}} \right)^{n}}\end{bmatrix}}},} & (1)\end{matrix}$ wherein [u_(x) ⁰ u_(y) ⁰]^(T) is value of the displacementu^(h) at the point P₀, and$a = {\left\lbrack {\frac{\partial u_{x}}{\partial x}❘_{P_{0}}{\frac{\partial u_{x}}{\partial y}❘_{P_{0}}{{\ldots\frac{\partial^{n}u_{x}}{\partial y^{n}}\frac{\partial u_{y}}{\partial x}}❘_{P_{0}}{\frac{\partial u_{y}}{\partial x}❘_{P_{0}}{\ldots\frac{\partial^{n}u_{y}}{\partial y^{n}}}}}}} \right\rbrack ❘_{P_{0}}}$is a displacement gradient from the first order to the nth order at thepoint P₀, where n is the order of Taylor expansion; step S12:determining support domain of each subdomain, and defining a weighteddiscrete L² norm J, wherein $\begin{matrix}{{J = {\left( {{Aa} + u_{0} - u_{m}} \right)^{T}{W\left( {{Aa} + u_{0} - u_{m}} \right)}}},} & (2)\end{matrix}$ $\begin{matrix}{A = \begin{bmatrix}{x_{1} - x_{0}} & {y_{1} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{1} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\{0} & 0 & {\ldots} & 0 & {x_{1} - x_{0}} & {y_{1} - y_{0}} & {\ldots} & {\frac{1}{n!}\left( {y_{1} - y_{0}} \right)^{n}} \\{x_{2} - x_{0}} & {y_{2} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{2} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & {x_{2} - x_{0}} & {y_{2} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{2} - y_{0}} \right)^{n}} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\{x_{m} - x_{0}} & {y_{m} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{m} - y_{0}} \right)^{n}} & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & {x_{m} - x_{0}} & {y_{m} - y_{0}} & \ldots & {\frac{1}{n!}\left( {y_{m} - y_{0}} \right)^{n}}\end{bmatrix}} & \end{matrix}$ $\begin{matrix}{{u_{0} = \left( \begin{bmatrix}{u_{x}^{0}} & u_{y}^{0} & u_{x}^{0} & u_{y}^{0} & \ldots & u_{x}^{0} & u_{y}^{0}\end{bmatrix}_{1 \times 2m} \right)^{T}},} & \end{matrix}$ $\begin{matrix}{{u_{m} = \begin{bmatrix}{u_{x}^{1}} & u_{y}^{1} & u_{x}^{2} & u_{y}^{2} & \ldots & u_{x}^{m} & u_{y}^{m}\end{bmatrix}^{T}},} & \end{matrix}$ $\begin{matrix}{{W = \begin{bmatrix}w_{x}^{1} & 0 & 0 & \ldots & 0 \\0 & w_{y}^{1} & \ddots & \ddots & \vdots \\0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & w_{x}^{m} & 0 \\0 & \ldots & 0 & 0 & w_{y}^{m}\end{bmatrix}},} & \end{matrix}$ (x_(i), y_(i)) is coordinates of a point P_(i), [u_(x)^(i) u_(y) ^(i)]^(T) is value of the displacement u^(h) at the pointP_(i), and [w_(x) ^(i) w_(y) ^(i)]^(T) is value of a weight function atthe point P_(i), (i=1, 2, 3, . . . m); step S13: expressing thedisplacement gradient a with point displacements according to theformula (2) and stationary conditions of the norm J, that is:$\begin{matrix}{{a = {Cu}_{E}},} & (3)\end{matrix}$ $\begin{matrix}{{{{wherein}C} = {\left( {A^{T}{WA}} \right)^{- 1}A^{T}{W\begin{bmatrix}I_{1} & I_{2}\end{bmatrix}}}},} & \end{matrix}$ $\begin{matrix}{{u_{E} = \begin{bmatrix}u_{x}^{0} & u_{y}^{0} & u_{x}^{1} & u_{y}^{1} & \ldots & u_{x}^{m} & u_{y}^{m}\end{bmatrix}^{T}},} & \end{matrix}$ $\begin{matrix}{{I_{1} = \begin{bmatrix}{- 1} & 0 \\0 & {- 1} \\{- 1} & 0 \\0 & {- 1} \\ \vdots & \vdots \\{- 1} & 0 \\0 & {- 1}\end{bmatrix}_{2m \times 2}},{{{{and}I_{2}} = \begin{bmatrix}1 & 0 & \ldots & 0 \\0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\0 & \ldots & 0 & 1\end{bmatrix}_{2m \times 2m}};}} & \end{matrix}$ step S14: obtaining the expression of the trial functionin each subdomain according to formulae (3) and (1); and u^(h) = Nu_(E)${ɛ^{h} = {\begin{bmatrix}ɛ_{11}^{h} \\ɛ_{22}^{h} \\{2ɛ_{12}^{h}}\end{bmatrix} = {{\begin{bmatrix}\frac{\partial}{\partial x_{1}} & 0 \\0 & \frac{\partial}{\partial x_{2}} \\\frac{\partial}{\partial x_{2}} & \frac{\partial}{\partial x_{1}}\end{bmatrix}u^{h}} = {Bu_{E}}}}},{\sigma^{h} = {\begin{bmatrix}\sigma_{11}^{h} \\\sigma_{22}^{h} \\\sigma_{12}^{h}\end{bmatrix} = {{Dɛ^{h}} = {DBu}_{E}}}}$ N is the shape function; D isthe stress-strain matrix; B is the derivative of shape function; ε^(h)is the strain and σ^(h) is the stress; step S15: obtaining theexpression of the test function by prescribing the test function in asame form as that of the trial function.
 16. The non-transitorycomputer-readable storage medium according to claim 15, whereindefinition of the support domain comprises any one of the followingmanners: defining the support domain of each point to contain allsubdomains having points adjacent to the corresponding point locatedtherein, or defining the support domain of each point to contain allsubdomains having points therein which are located in a circle of aparticular size and with the corresponding point at the center; and if acrack exists between two neighboring subdomains, the point located onone side of the crack is not contained in the support domain of thepoint on the other side of the crack.
 17. The non-transitorycomputer-readable storage medium according to claim 13, wherein thenumerical flux corrections comprises interior penalty numerical fluxcorrections.